{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/NeuromatchAcademy/course-content/blob/master/tutorials/W1D3_ModelFitting/W1D3_Tutorial5.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "# Neuromatch Academy: Week 1, Day 3, Tutorial 5\n",
    "# Model Selection: Bias-variance trade-off\n",
    "\n",
    "**Content creators**: Pierre-Étienne Fiquet, Anqi Wu, Alex Hyafil with help from Ella Batty\n",
    "\n",
    "**Content reviewers**: Lina Teichmann, Patrick Mineault, Michael Waskom\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "#Tutorial Objectives\n",
    "\n",
    "This is Tutorial 5 of a series on fitting models to data. We start with simple linear regression, using least squares optimization (Tutorial 1) and Maximum Likelihood Estimation (Tutorial 2). We will use bootstrapping to build confidence intervals around the inferred linear model parameters (Tutorial 3). We'll finish our exploration of regression models by generalizing to multiple linear regression and polynomial regression (Tutorial 4). We end by learning how to choose between these various models. We discuss the bias-variance trade-off (Tutorial 5) and Cross Validation for model selection (Tutorial 6).\n",
    "\n",
    "In this tutorial, we will learn about the bias-variance tradeoff and see it in action using polynomial regression models.\n",
    "\n",
    "Tutorial objectives:\n",
    "\n",
    "* Understand difference between test and train data\n",
    "* Compare train and test error for models of varying complexity\n",
    "* Understand how bias-variance tradeoff relates to what model we choose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 517
    },
    "colab_type": "code",
    "outputId": "8d20c293-126c-47d3-bdce-3c108f1da640"
   },
   "outputs": [],
   "source": [
    "#@title Video 1: Bias Variance Tradeoff\n",
    "from IPython.display import YouTubeVideo\n",
    "video = YouTubeVideo(id=\"NcUH_seBcVw\", width=854, height=480, fs=1)\n",
    "print(\"Video available at https://youtube.com/watch?v=\" + video.id)\n",
    "video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "#@title Figure Settings\n",
    "%config InlineBackend.figure_format = 'retina'\n",
    "plt.style.use(\"https://raw.githubusercontent.com/NeuromatchAcademy/course-content/master/nma.mplstyle\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "#@title Helper functions\n",
    "def ordinary_least_squares(x, y):\n",
    "  \"\"\"Ordinary least squares estimator for linear regression.\n",
    "\n",
    "  Args:\n",
    "    x (ndarray): design matrix of shape (n_samples, n_regressors)\n",
    "    y (ndarray): vector of measurements of shape (n_samples)\n",
    "\n",
    "  Returns:\n",
    "    ndarray: estimated parameter values of shape (n_regressors)\n",
    "  \"\"\"\n",
    "\n",
    "  return np.linalg.inv(x.T @ x) @ x.T @ y\n",
    "\n",
    "def make_design_matrix(x, order):\n",
    "  \"\"\"Create the design matrix of inputs for use in polynomial regression\n",
    "\n",
    "  Args:\n",
    "    x (ndarray): input vector of shape (n_samples)\n",
    "    order (scalar): polynomial regression order\n",
    "\n",
    "  Returns:\n",
    "    ndarray: design matrix for polynomial regression of shape (samples, order+1)\n",
    "  \"\"\"\n",
    "\n",
    "  # Broadcast to shape (n x 1) so dimensions work\n",
    "  if x.ndim == 1:\n",
    "    x = x[:, None]\n",
    "\n",
    "  #if x has more than one feature, we don't want multiple columns of ones so we assign\n",
    "  # x^0 here\n",
    "  design_matrix = np.ones((x.shape[0],1))\n",
    "\n",
    "  # Loop through rest of degrees and stack columns\n",
    "  for degree in range(1, order+1):\n",
    "      design_matrix = np.hstack((design_matrix, x**degree))\n",
    "\n",
    "  return design_matrix\n",
    "\n",
    "\n",
    "def solve_poly_reg(x, y, max_order):\n",
    "  \"\"\"Fit a polynomial regression model for each order 0 through max_order.\n",
    "\n",
    "  Args:\n",
    "    x (ndarray): input vector of shape (n_samples)\n",
    "    y (ndarray): vector of measurements of shape (n_samples)\n",
    "    max_order (scalar): max order for polynomial fits\n",
    "\n",
    "  Returns:\n",
    "    dict: fitted weights for each polynomial model (dict key is order)\n",
    "  \"\"\"\n",
    "\n",
    "  # Create a dictionary with polynomial order as keys, and np array of theta\n",
    "  # (weights) as the values\n",
    "  theta_hats = {}\n",
    "\n",
    "  # Loop over polynomial orders from 0 through max_order\n",
    "  for order in range(max_order+1):\n",
    "\n",
    "    X = make_design_matrix(x, order)\n",
    "    this_theta = ordinary_least_squares(X, y)\n",
    "\n",
    "    theta_hats[order] = this_theta\n",
    "\n",
    "  return theta_hats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 1: Train vs test data\n",
    "\n",
    " The data used for the fitting procedure for a given model is the **training data**. In tutorial 4, we computed MSE on the training data of our polynomial regression models and compared training MSE across models. An additional important type of data is **test data**. This is held-out data that is not used (in any way) during the fitting procedure. When fitting models, we often want to consider both the train error (the quality of prediction on the training data) and the test error (the quality of prediction on the test data) as we will see in the next section.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "We will generate some noisy data for use in this tutorial using a similar process as in Tutorial 4.However, now we will also generate test data. We want to see how our model generalizes beyond the range of values see in the training phase. To accomplish this, we will generate x from a wider range of values ([-3, 3]). We then plot the train and test data together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "colab_type": "code",
    "outputId": "2d24db63-27eb-46a1-9182-b21b64d622bd"
   },
   "outputs": [],
   "source": [
    "#@title\n",
    "#@markdown Execute this cell to simulate both training and test data\n",
    "\n",
    "### Generate training data\n",
    "np.random.seed(0)\n",
    "n_train_samples = 50\n",
    "x_train = np.random.uniform(-2, 2.5, n_train_samples) # sample from a uniform distribution over [-2, 2.5)\n",
    "noise = np.random.randn(n_train_samples) # sample from a standard normal distribution\n",
    "y_train =  x_train**2 - x_train - 2 + noise\n",
    "\n",
    "### Generate testing data\n",
    "n_test_samples = 20\n",
    "x_test = np.random.uniform(-3, 3, n_test_samples) # sample from a uniform distribution over [-2, 2.5)\n",
    "noise = np.random.randn(n_test_samples) # sample from a standard normal distribution\n",
    "y_test =  x_test**2 - x_test - 2 + noise\n",
    "\n",
    "## Plot both train and test data\n",
    "fig, ax = plt.subplots()\n",
    "plt.title('Training & Test Data')\n",
    "plt.plot(x_train, y_train, '.', markersize=15, label='Training')\n",
    "plt.plot(x_test, y_test, 'g+', markersize=15, label='Test')\n",
    "plt.legend()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('y');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Section 2: Bias-variance tradeoff\n",
    "\n",
    "Finding a good model can be difficult. One of the most important concepts to keep in mind when modeling is the **bias-variance tradeoff**. \n",
    "\n",
    "**Bias** is the difference between the prediction of the model and the corresponding true output variables you are trying to predict. Models with high bias will not fit the training data well since the predictions are quite different from the true data. These high bias models are overly simplified - they do not have enough parameters and complexity to accurately capture the patterns in the data and are thus **underfitting**.\n",
    "\n",
    "\n",
    "**Variance** refers to the variability of model predictions for a given input. Essentially, do the model predictions change a lot with changes in the exact training data used? Models with high variance are highly dependent on the exact training data used - they will not generalize well to test data. These high variance models are **overfitting** to the data.\n",
    "\n",
    "In essence:\n",
    "\n",
    "* High bias, low variance models have high train and test error.\n",
    "* Low bias, high variance models have low train error, high test error\n",
    "* Low bias, low variance models have low train and test error\n",
    "\n",
    "\n",
    "As we can see from this list, we ideally want low bias and low variance models! These goals can be in conflict though - models with enough complexity to have low bias also tend to overfit and depend on the training data more. We need to decide on the correct tradeoff.\n",
    "\n",
    "In this section, we will see the bias-variance tradeoff in action with polynomial regression models of different orders.\n",
    "\n",
    "Graphical illustration of bias and variance.\n",
    "(Source: http://scott.fortmann-roe.com/docs/BiasVariance.html)\n",
    "\n",
    "![bias-variance](https://www.cs.cornell.edu/courses/cs4780/2018fa/lectures/images/bias_variance/bullseye.png) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "We will first fit polynomial regression models of orders 0-5 on our simulated training data just as we did in Tutorial 4. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "#@title\n",
    "#@markdown Execute this cell to estimate theta_hats\n",
    "max_order = 5\n",
    "theta_hats = solve_poly_reg(x_train, y_train, max_order)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "### Exercise 1: Compute and compare train vs test error\n",
    "\n",
    "We will use MSE as our error metric again. Compute MSE on training data ($x_{train},y_{train}$) and test data ($x_{test}, y_{test}$) for each polynomial regression model (orders 0-5). Since you already developed code in T4 Exercise 4 for evaluating fit polynomials, we have ported that here into the function ``evaluate_poly_reg`` for your use.\n",
    "\n",
    "*Please think about after completing exercise before reading the following text! Do you think the order 0 model has high or low bias? High or low variance? How about the order 5 model?*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code"
   },
   "outputs": [],
   "source": [
    "def evaluate_poly_reg(x, y, theta_hats, max_order):\n",
    "    \"\"\" Evaluates MSE of polynomial regression models on data\n",
    "\n",
    "    Args:\n",
    "      x (ndarray): input vector of shape (n_samples)\n",
    "      y (ndarray): vector of measurements of shape (n_samples)\n",
    "      theta_hats (dict):  fitted weights for each polynomial model (dict key is order)\n",
    "      max_order (scalar): max order of polynomial fit\n",
    "\n",
    "    Returns\n",
    "      (ndarray): mean squared error for each order, shape (max_order)\n",
    "    \"\"\"\n",
    "\n",
    "    mse = np.zeros((max_order + 1))\n",
    "    for order in range(0, max_order + 1):\n",
    "      X_design = make_design_matrix(x, order)\n",
    "      y_hat = np.dot(X_design, theta_hats[order])\n",
    "      residuals = y - y_hat\n",
    "      mse[order] = np.mean(residuals ** 2)\n",
    "\n",
    "    return mse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 447
    },
    "colab_type": "code",
    "outputId": "c242cb62-5e84-42a8-d01d-14ad4ce6475c"
   },
   "outputs": [],
   "source": [
    "def compute_mse(x_train,x_test,y_train,y_test,theta_hats,max_order):\n",
    "  \"\"\"Compute MSE on training data and test data.\n",
    "\n",
    "  Args:\n",
    "    x_train(ndarray): training data input vector of shape (n_samples)\n",
    "    x_test(ndarray): test data input vector of shape (n_samples)\n",
    "    y_train(ndarray): training vector of measurements of shape (n_samples)\n",
    "    y_test(ndarray): test vector of measurements of shape (n_samples)\n",
    "    theta_hats(dict): fitted weights for each polynomial model (dict key is order)\n",
    "    max_order (scalar): max order of polynomial fit\n",
    "\n",
    "  Returns:\n",
    "    ndarray, ndarray: MSE error on training data and test data for each order\n",
    "  \"\"\"\n",
    "\n",
    "  #######################################################\n",
    "  ## TODO for students: calculate mse error for both sets\n",
    "  ## Hint: look back at tutorial 5 where we calculated MSE\n",
    "  # Fill out function and remove\n",
    "  raise NotImplementedError(\"Student excercise: calculate mse for train and test set\")\n",
    "  #######################################################\n",
    "\n",
    "  mse_train = ...\n",
    "  mse_test = ...\n",
    "\n",
    "  return mse_train, mse_test\n",
    "\n",
    "\n",
    "#Uncomment below to test your function\n",
    "# mse_train, mse_test = compute_mse(x_train, x_test, y_train, y_test, theta_hats, max_order)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "width = .35\n",
    "\n",
    "# ax.bar(np.arange(max_order + 1) - width / 2, mse_train, width, label=\"train MSE\")\n",
    "# ax.bar(np.arange(max_order + 1) + width / 2, mse_test , width, label=\"test MSE\")\n",
    "\n",
    "ax.legend()\n",
    "ax.set(xlabel='Polynomial order', ylabel='MSE', title ='Comparing polynomial fits');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "both",
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 430
    },
    "colab_type": "code",
    "outputId": "959f79ef-68d2-474d-9e3e-0666a9f8e83e"
   },
   "outputs": [],
   "source": [
    "# to_remove solution\n",
    "\n",
    "def compute_mse(x_train,x_test,y_train,y_test,theta_hats,max_order):\n",
    "  \"\"\"Compute MSE on training data and test data.\n",
    "\n",
    "  Args:\n",
    "    x_train(ndarray): training data input vector of shape (n_samples)\n",
    "    x_test(ndarray): test vector of shape (n_samples)\n",
    "    y_train(ndarray): training vector of measurements of shape (n_samples)\n",
    "    y_test(ndarray): test vector of measurements of shape (n_samples)\n",
    "    theta_hats(dict): fitted weights for each polynomial model (dict key is order)\n",
    "    max_order (scalar): max order of polynomial fit\n",
    "\n",
    "  Returns:\n",
    "    ndarray, ndarray: MSE error on training data and test data for each order\n",
    "  \"\"\"\n",
    "\n",
    "  mse_train = evaluate_poly_reg(x_train, y_train, theta_hats, max_order)\n",
    "  mse_test = evaluate_poly_reg(x_test, y_test, theta_hats, max_order)\n",
    "\n",
    "  return mse_train, mse_test\n",
    "\n",
    "\n",
    "mse_train, mse_test = compute_mse(x_train, x_test, y_train, y_test, theta_hats, max_order)\n",
    "\n",
    "with plt.xkcd():\n",
    "  fig, ax = plt.subplots()\n",
    "  width = .35\n",
    "\n",
    "  ax.bar(np.arange(max_order + 1) - width / 2, mse_train, width, label=\"train MSE\")\n",
    "  ax.bar(np.arange(max_order + 1) + width / 2, mse_test , width, label=\"test MSE\")\n",
    "\n",
    "  ax.legend()\n",
    "  ax.set(xlabel='Polynomial order', ylabel='MSE', title ='Comparing polynomial fits');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "As we can see from the plot above, more complex models (higher order polynomials) have lower MSE for training data. The overly simplified models (orders 0 and 1) have high MSE on the training data. As we add complexity to the model, we go from high bias to low bias. \n",
    "\n",
    "The MSE on test data follows a different pattern. The best test MSE is for an order 2 model - this makes sense as the data was generated with an order 2 model. Both simpler models and more complex models have higher test MSE. \n",
    "\n",
    "So to recap:\n",
    "\n",
    "Order 0 model: High bias, low variance\n",
    "\n",
    "Order 5 model: Low bias, high variance\n",
    "\n",
    "Order 2 model: Just right, low bias, low variance\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Summary\n",
    "\n",
    "- Training data is the data used for fitting, test data is held-out data.\n",
    "- We need to strike the right balance between bias and variance. Ideally we want to find a model with optimal model complexity that has both low bias and low variance\n",
    " - Too complex models have low bias and high variance.\n",
    " - Too simple models have high bias and low variance."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "**Note**\n",
    " - Bias and variance are very important concepts in modern machine learning, but it has recently been observed that they do not necessaruly trade off (see for example the phenomenon and theory of \"double descent\")\n",
    "\n",
    "**Further readings:**\n",
    "- [The elements of statistical learning](https://web.stanford.edu/~hastie/ElemStatLearn/) by Hastie, Tibshirani and Friedman"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "---\n",
    "# Appendix\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text"
   },
   "source": [
    "## Bonus Exercise\n",
    "\n",
    "Prove the bias-variance decomposition for MSE\n",
    "\n",
    "$$\n",
    "\\mathrm{E}_{x}\\left[(y-\\hat{y}(x ; \\theta))^{2}\\right]=\\left(\\operatorname{Bias}_{x}[\\hat{y}(x ; \\theta)]\\right)^{2}+\\operatorname{Var}_{x}[\\hat{y}(x ; \\theta)]+\\sigma^{2}\n",
    "$$where\n",
    "$$\\operatorname{Bias}_{x}[\\hat{y}(x ; \\theta)]=\\mathrm{E}_{x}[\\hat{y}(x ; \\theta)]-y\n",
    "$$and\n",
    "$$\\operatorname{Var}_{x}[\\hat{y}(x ; \\theta)]=\\mathrm{E}_{x}\\left[\\hat{y}(x ; \\theta)^{2}\\right]-\\mathrm{E}_{x}[\\hat{y}(x ; \\theta)]^{2}\n",
    "$$\n",
    "\n",
    "Hint: use $$\\operatorname{Var}[X]=\\mathrm{E}\\left[X^{2}\\right]-(\\mathrm{E}[X])^{2}$$"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "include_colab_link": true,
   "name": "W1D3_Tutorial5",
   "provenance": [],
   "toc_visible": true
  },
  "kernel": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  },
  "toc-autonumbering": true
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
